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Drug gradients are believed to play an important role in the evolution of drug resistance, from 
antibiotics to cancer. We use a statistical physics model to study the evolution of a population of 
malignant cells exposed to drug gradients, where drug resistance emerges via multiple mutations. 
We show that a non-uniform drug concentration can strongly accelerate the emergence of resistance 
when the mutational pathway involves a long sequence of mutants with increasing resistance, but 
slows it down if the pathway is short or crosses a fitness valley. These predictions can be verified 
experimentally, and have the potential to improve strategies to combat the emergence of resistance. 

PACS numbers: 02.50.Ey, 05.70.Fh, 05.70.Ln, 64.60.-i 



The evolution of drug resistance is an urgent problem 
in the treatment of diseases, from bacterial infections to 
cancer. Attempts to address this problem include the 
characterization of mutational pathways leading to resis- 
tance 2J, as well as theoretical [3-8J and experimen- 
tal [OHTlJ studies of the emergence of resistance under 
different treatment regimens. These studies often as- 
sume a spatially uniform drug concentration. However, 
in many clinical situations drug concentrations vary in 
space [T2I |13|, for example where malignant cells form 
less drug-permeable layers such as bacterial biofilms [MJ 
or tumour stromas [15j. Recent experimental work [16] 
suggests that the evolution of antibiotic resistance in bac- 
terial populations can be greatly accelerated if the antibi- 
otic concentration is spatially non-uniform. 

It is often observed that several mutations are required 
to obtain maximal resistance to a drug [T1[5J[TB]. In some 
cases, fitness (i.e. drug resistance) increases steadily 
along the mutational pathway to full resistance [l]; in 
other cases, epistatic interactions between mutations may 
result in less fit intermediate genotypes (fitness "valleys") 
fTTtilQ] . The role of mutational pathways in controlling 
evolutionary dynamics has been studied in the quasis- 
pecies model [20j [21] and in models of cancer progression 
[551 [23j. That work, however, does not take into account 
the effects of spatial structure. On the other hand, in 
models without complex evolutionary pathways, it is well 
known that spatial structure can increase genetic diver- 
sity, the rate of evolutionary diversification [24H28], and 
the rate of viral drug resistance [29j; indeed, in a broader 
statistical physics context, spatial structure plays a key 
role in many theoretical studies of evolving populations 

[21I3DH35]. 

Here we use a statistical physics model to show that, 
in the presence of drug gradients, a population evolves 
drug resistance in a sequence of waves of increasingly 
better adapted mutants that extend its range in a step- 
wise manner. In contrast, for a uniformly distributed 
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Figure 1: Simulation snapshots for the cases of a uniform 
and an exponentially increasing drug concentration (left and 
right, respectively). Blue lines show the drug concentration 
(left axes), while the colors represent the populations of the 
different genotypes (right axes). Parameter values are K = 
100, L = 500, M = 6,/i = 5 X 10"^ c™^^ = 4^-\ and the 
drug concentration c = 0.3 (left panel) and c(x) = e^^ — 1 
with a — 0.06 (right panel). For corresponding movies see 



drug, resistant mutants evolve at random positions and 
spread over the entire environment. The rate of evolu- 
tion depends crucially on the mutational pathway lead- 
ing to drug resistance. If tolerance to the drug increases 
monotonously along the pathway, drug gradients can sig- 
nificantly accelerate the evolution by increasing selection 
at the population's edge. However, if the pathway crosses 
a fitness valley, evolution of resistance may actually be 
slowed down by a non-uniform drug concentration, as a 
result of a reduced rate of "stochastic tunnelling" due to a 
smaller population size. We also suggest how our predic- 
tions could be verified experimentally and used to infer 
the structure of the mutational pathway to resistance. 

The model. We consider a growing population of cells 
which mutate between M possible genotypes, with differ- 
ent levels of resistance to a drug. To model the effects of 
spatial heterogeneity, the population is assumed to reside 
within a chain of L connected microhabitats, which may 
contain different concentrations of the drug. These dis- 
crete microhabitats might represent connected chambers 
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in a microfiuidic experiment fTB]; in the limit of small 
microhabitats, our model represents a population grow- 
ing in continuous space. Within a given microhabitat i 
the population is assumed to be well-mixed, with a fixed 
carrying capacity K; cells of genotype m replicate at rate 
(j)m{ci){l—Ni/ K) where Ni is the total population of cells 
in microhabitat q is the drug concentration and (j)m{ci) 
is the growth rate of genotype m, which depends on the 
local drug concentration. Upon replication, cells mutate 
with probability we shall mainly consider the case of 
an unbranched mutational pathway to drug resistance, 
such that genotype m mutates only into genotypes mil 
(without any bias). Cells migrate between microhabitats 
i and z ± 1 at rate h/2 and die at a fixed rate d. 

A key feature of our model is the fact that differ- 
ent genotypes show different levels of drug resistance: 
genotype 1 is least resistant while genotype M is most 
resistant. The minimal inhibitory concentration (MIC) 
denotes the drug concentration at which genotype m 
ceases to be able to grow; this is embodied in our model 

in the assumption that ^rn(c) = max |o, 1 — (c/c™^*^)^|. 
This choice is inspired by Ref. [37J (see also [38j). 

We study this model using kinetic Monte Carlo sim- 
ulations [39j which are initiated with Ni = K cells of 
genotype 1 in microhabitat 1 and with all the other mi- 
crohabitats empty, so that the population colonizes the 
space during the simulation. We define the units of time 
in our simulations by fixing the maximal growth rate 
^m(O) = 1 and the units of drug concentration by fix- 
ing 1, and we set b 0.1, d 0.1, and K 100 
(see [38j for a detailed discussion of the choice of the 
parameters). We use values of /i and L such that the 
number of mutants per generation emerging anywhere in 
the environment is typically small, jj^KL < 1 (see fig- 
ure captions for details). Our results are thus relevant 
for not-too-large populations, for which stochastic effects 
play an important role. 

Monotonically increasing MIC. We first consider a 
pathway to resistance for which the M genotypes have 
increasing levels of drug resistance - i.e. c^'^ > c^^^^ for 
all m > 1, as depicted in Fig. [2]:. In particular, we set 
M = 6 and = 4^-^ the ratio cf'^/cf'^ ^ 10^ be- 
tween fully resistant and wild-type cells is consistent with 
experimentally determined values [Ij (see also [38j). We 
compare the "homogeneous" case where the drug concen- 
tration is uniform (q = c) with the "spatially heteroge- 
neous" case where the drug concentration q = exp(ai) — 1 
is non-uniform and increases exponentially from left to 
right with steepness a (see solid lines in Fig. [T]) . 

The emergence of drug resistance occurs very differ- 
ently in the homogeneous and heterogeneous cases, as 
illustrated in the simulation snapshots of Fig. [T] In the 
homogeneous case (Fig. [l] left), genotype 1 (blue) first 
spreads to fill the entire space, then mutants of geno- 
type 2 (green) emerge at random locations; these spread 
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Figure 2: Average time to resistance f for homogeneous 
((a,d), red circles) and heterogeneous ((b,e), blue circles) drug 
concentrations, for M = 6,L = 500, = 100, and /x = 5 x 
10~^. Panels (a,d): exponentially increasing MIC (c), panels 
(b,e): fitness valley (f). For the heterogeneous case (b,e), 
the red dashed lines show the minimal value of f obtained 
for the homogeneous case (i.e. the minimum of t(c) from 
(a,d)) . The black lines show the theoretical predictions of Eqs. 
Q (exponential MIC) and Eqs. ^ (valley MIC). The 

insets show simulation shapshots taken just before the first 
occurence of genotype m = 6, for two values of a (indicated 
by arrows). 



because they have a growth advantage in the drug en- 
vironment, before giving rise to genotype 3 (red), etc. 
[47J. In contrast, in the heterogeneous case, waves of 
increasingly better-adapted mutants invade the environ- 
ment from left to right in a step- wise manner. Each geno- 
type colonizes the space only up to a well-defined spatial 
boundary, where the local drug concentration approaches 
its MIC and the population forms a stationary "front". 
Better-adapted mutants then emerge at the edge of this 
front to further colonize the space. 

To quantify the effect of a spatial drug distribution on 
the emergence of drug resistance, we plot in Fig. |2^,b the 
mean time f to emergence of full drug resistance - i.e. 
the time to emergence of a mutant with m = M = 6, 
averaged over surviving populations. For a homoge- 
neous drug concentration q = c, f decreases as the drug 
concentration c increases (Fig. |2|i). Resistance emerges 
fastest when c approaches the MIC of genotype 1 (de- 
fined to be unity). For drug concentrations above C]"^^, 
however, the population does not evolve resistance be- 
cause the initial genotype cannot reproduce. For the 
non-uniform drug concentration (Fig. [2]3), f varies non- 
monotonically as a function of the steepness a, with a 
minimum at a ~ 0.01. This minimum arises because for 
very small a, little drug is present, reducing the selection 
pressure for the evolution of resistant mutants, whereas 
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for very large a, the fronts become narrow, reducing the 
size of the "zone" in which new resistant mutants can 
emerge (see snapshots in Fig. |2]3). 

Importantly, for almost all values of a and c, resis- 
tance emerges faster for the non-uniform drug concen- 
tration than for the uniform case; the minimal value of 
r(a) in the non-uniform case is smaller by an order of 
magnitude than the minimal value of r(c) in the uni- 
form case (dashed line in Fig. [2]3). Thus, if the MIC in- 
creases monotonically along the pathway to resistance, a 
non-uniform drug concentration carries the potential for 
much faster evolution of drug resistance than is possible 
in a uniform drug concentration. 

These results can be rationalized using simple physical 
arguments. For the non-uniform drug concentration we 
consider separately the colonization of space by a subpop- 
ulation of cells of genotype m and the stochastic emer- 
gence of mutants of genotype m + 1 within this subpop- 
ulation. In the continuous approximation (valid for large 
L and a <C 1) and assuming a unit distance between the 
habitats, the expansion of a wave of mutants of genotype 
m is described by the Fisher-KPP equation [40j: 
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where x = (j)rn = ^m(c(^)) and = Nj^^x^t) de- 
notes the population of genotype m. If c{x) <C c] 
(prn ^ 1 and Eq. ([T]) describes a Fisher wave prop- 
agating with speed v ^ y^26(1 — d) [40j. The wave 
stops when it reaches the point where c{x) ^ c^^^; 
for small b the stationary solution of Eq. ([T]) reads 
N^{x) = K[l — d/(j)rn{c{x))]^ which decays to zero at 
= {l/a)\n{c^'''^/(l^dj ^ 1). Assuming that the 
wave of mutants of genotype m emerges at x^_i (i.e. 
at the stationary front of the preceding wave), the time 
it takes to reach its stationary state is then T^^^® ~ 
{x*^ - x*^_,)/v, with Tr"' « xl/v. 

Once the stationary population of genotype m is estab- 
lished, the waiting time before a new wave of mutants of 
genotype m + 1 arises can be expressed for low mutation 
rates as the inverse of the total rate at which mutants 
establish in the population. 
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Here /i/2 is the probability to mutate from genotype m to 
m + 1, r{x) ^ d is the rate of reproduction in the steady 
state and Pfix = (0m+i — ^m)/^m+i is the probability of 
fixation of genotype m + 1; this is a standard result [41j . 
The mean total time until the first cell of genotype M 




Figure 3: Average time r to full resistance (left) and its co- 
efficient of variation Cv — yr^^—f^/f (right) as a function of 
the mutational pathway length M for homogeneous (c = 0.9; 
triangles) and heterogeneous {a = 0.07; pluses) drug concen- 
tration. In both cases L = 300, K = 100, and /j, = 10~^. Blue 
lines are theoretical predictions for the heterogeneous case: 
the mean f is calculated numerically from Eqs. ^ and ([3| 
and Cv = A^/M / {M + B) with A and B fitted to data. 
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Note that we neglect the short time for the first cell of 
genotype M to occur after genotype M — 1 has reached 
its stationary value, and we assume strong selection, so 
that backward mutations do not fix. Equation ([3| de- 
composes the time to resistance into the independent 
contributions of each wave of mutants. The value of 
f calculated from Eqs. ^ and ([3| is in good quanti- 
tive agreement with our simulation results (black line 
in Fig. Since each successive emergence of a new 

genotype m can be modelled as a Poisson process, we 
expect that in a given experiment, the measured value of 
r will be equal to a constant contribution ^^Z^ ^m^^^ 
plus the sum of M — 2 exponentially-distributed Pois- 
son waiting times with means T^^^. For our choice of 
MICs and drug distribution, T^""^ and T^^^^ are approx- 
imately independent of m for large m, because then the 
stationary waves N^{x) just differ by a constant spatial 
shift [42j. Hence for our model we expect r to follow ap- 
proximately a shifted Erlang distribution with the mean 
f growing linearly with the length of the pathway M, 
and its coefficient of variation Cj, = \/r^ — f^/f scaling 
as >/M / (const + M) . Figure [s] shows that this is in- 
deed the case in our simulations, for pathways of length 
M = 2, . . . , 20. This prediction breaks down, however, 
for M > merit ^ aL/ln(4) + 2 ^ 17, because the pop- 
ulation hits the right boundary before genotype M has 
evolved [48j. 

A different argument applies in the case where the 
drug is uniformly distributed, c{x) = c. For short path- 
ways (small M), resistance evolves faster in the uniform 
environment than in the non-uniform one (Fig. |3| left 
panel) , since the population size from which mutants can 
emerge is larger. However, beyond a critical pathway 
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length (M > 4 in our simulations), resistance emerges 
more slowly in the uniform environment, even though 
the population size is larger. In the uniform environ- 
ment, new genotypes must establish in an already fully 
colonized space at fixed drug concentration. The time to 
fixation of genotype m + 1 is expected to scale with the 
fitness difference as — ^^)~^, where 7 > depends 

on whether mutations are rare or frequent f43]. In our 
model, the fitness advantage of genotype m + 1 over m is 
0m+i - </>m = c2((l/c-i-)2 - which dccrcascs 

for successive genotypes, because c^^^ increases with m. 
Therefore, successive genotypes take longer to fix, lead- 
ing to a super-linear increase of f(M) with M. This 
result is also valid for other functional forms of ^m(c), as 
long as resistance increases with m. For long pathways, 
the scaling changes again: here, selection for successive 
mutants becomes so weak that backward mutations can 
establish and neutral drift dominates. The evolution can 
then be approximated by an unbiased random walk in 
genotype space; f is then given by the mean first passage 
time ^ M of a walker constrained to positive m, which 
scales as ~ [44J, with coefficient of variation Cy inde- 
pendent of M (see Fig.jsj right panel). Note that here we 
have assumed that the MICs of intermediate genotypes 
remain fixed as M changes; one can also show [45j that 
the same general conclusions hold (including the hetero- 
geneous case) if c^^^ is scaled with M so as to keep the 
MIC of the most resistant genotype constant. 

Fitness valley. We now contrast these results with 
the case where the pathway to resistance passes through 
a "fitness valley" - i.e. one of the intermediate geno- 
types m has a lower MIC (is less drug-resistant) than its 
neighbouring genotypes m — 1 and m + 1. This scenario 
can arise due to epistatic interactions between mutations 
[l71[18j, such that two mutations are required to gain a 
particular fitness benefit. In particular, in our model we 
set Cg^^^ = 3.5, keeping all the other c^^^ = 4^~^ as be- 
fore, so that cf'"" > cf''' < cf'"", as depicted in Fig. 
Figure [2]i,e shows that the presence of the fitness valley 
has a dramatic effect on the time to resistance in the het- 
erogeneous environment: f{a) now rises steeply with a. 
Crucially, the shortest time to resistance in the heteroge- 
neous environment is now comparable to that in the ho- 
mogeneous environment, mma{T{a)) ~ minc(r(c)), and 
r{a) > minc(r(c)) for almost all values of a. Thus, when 
the pathway to resistance contains a fitness valley, a non- 
uniform drug concentration does not speed up, and may 
well slow down the emergence of resistance. 

To understand this effect, we argue that the rate- 
limiting step in the evolutionary process is the "tun- 
nelling" through the fitness valley [46j: mutants of geno- 
type 4 arise from the population of genotype 2 via short- 
lived mutants of genotype 3 which do not reach fixation. 
The tunnelling rate has been calculated for well- 
mixed populations in Ref. [46j (Eq. (2) therein). Apply- 
ing this result to the case of uniform drug distribution we 



obtain f-i ^ rA^2(/i/2)^(Pfix/5) where N2 ^ LK{1 - d) 
is the population size of genotype 2, s — (^2 — ^3)/ ^2 
is the selective advantage of genotype 2 over genotype 3, 
^fix = (04 — 02)/02 is the fixation probability of geno- 
type 4 and r is the growth rate which in the steady state 
equals the death rate d. For our choice of (\){c) and d 
this gives 
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f ^ 1.23/(ViV2) = \.2Zl{d^^LK(\ - d)), (4) 

which is independent of c. Equation Q is in good agree- 
ment with our simulation results (black line in Fig. [2]i). 
Extending this approach to the heterogeneous case, we 
integrate over the steady-state population density of 
genotype 2: 
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This result agrees well with our simulation results for 
the non- uniform drug concentration (Fig. [2^). The in- 
crease in f with the steepness of the drug concentration 
profile a occurs because the domain occupied by popu- 
lation 2 decreases as a increases; the non-uniform drug 
concentration decreases the steady state population size 
of genotype 2, reducing the pool of cells from which mu- 
tants of genotype 4 can emerge and slowing down the 
evolution of resistance. 

Conclusion. Our results show that the mutational 
pathway to drug resistance plays a crucial role in de- 
termining the effect of a spatial drug distribution on the 
time to evolve resistance. If fitness (i.e. level of drug 
resistance) increases monotonically along the mutational 
pathway, a non-uniform drug concentration can greatly 
accelerate the evolution of resistance, by a factor that 
increases dramatically with the length of the pathway. 
However, for short pathways, or those involving a fitness 
valley, a non-uniform drug concentration does not speed 
up the evolution of resistance - indeed, it may actually 
slow it down. 

Our predictions can be verified experimentally. Recent 
microfluidic experiments have shown that gradients of 
the antibiotic ciprofloxacin greatly accelerate the emer- 
gence of resistance of the bacterium E. coli [16j. Al- 
though the mutational pathway in this case is not known, 
our results suggest that it is likely to be monotonic. Fur- 
thermore, we predict that repeating the experiments us- 
ing cefotaxime (monotonic pathway [Ij) should produce 
similar results, but that for streptomycin, which has a 
fitness valley [T7l [18], resistance should emerge faster in 
a uniform drug concentration. This procedure could be 
generalized to infer the characteristics of unknown muta- 
tional pathways, by comparing the times to resistance of 
cells in a microfluidic device, for different drug concen- 
trations and gradients. 
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